!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==============================================================================================================================================|
subroutine flux (ro, mx, my, mz, bx, by, bz, en, fro, fmx, fmy, fmz, fbx, fby, fbz, fen, nx, gm)
!==============================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx

   real(8), intent(in) :: gm
   real(8), dimension(nx), intent(in) :: ro, mx, my, mz, bx, by, bz, en
   real(8), dimension(nx), intent(inout) :: fro, fmx, fmy, fmz, fbx, fby, fbz, fen

   integer(4) :: i

   real(8) :: gmm1, vx, vy, vz, b2, v2, pr

!----------------------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0

   do i = 1, nx
      vx = mx(i) / ro(i)
      vy = my(i) / ro(i)
      vz = mz(i) / ro(i)
      b2 = bx(i) * bx(i) + by(i) * by(i) + bz(i) * bz(i)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en(i) - v2 * ro(i) / 2.0d0 - b2 / 2.0d0) * gmm1
      fro(i) = mx(i)
	   fmx(i) = vx * mx(i) + pr + b2 / 2.0d0 - bx(i) * bx(i)
	   fmy(i) = vx * my(i) - bx(i) * by(i)
	   fmz(i) = vx * mz(i) - bx(i) * bz(i)
	   fbx(i) = 0.0d0
	   fby(i) = vx * by(i) - bx(i) * vy
	   fbz(i) = vx * bz(i) - bx(i) * vz
	   fen(i) = (en(i) + pr + b2 / 2.0d0) * vx - (bx(i) * vx + by(i) * vy + bz(i) * vz) * bx(i)
   enddo

!----------------------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine flux
